Approaching condensed matter ground states from below 
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We present a generic approach to the condensed matter ground state problem which is complemen- 
tary to variational techniques and works directly in the thermodynamic limit. Relaxing the ground 
state problem, we obtain semidefinite programs (SDP). These can be solved efficiently, yielding 
lower bounds to the ground state energy and approximations to the few-particle Green's functions. 
As the method is applicable for all particle statistics, it represents in particular a novel route for 
the study of strongly correlated fermionic and frustrated spin systems in D > 1 spatial dimensions. 
It is demonstrated for the XXZ model and the Hubbard model of spinless fermions. The results 
are compared against exact solutions, Quantum Monte Carlo, and Anderson bounds, showing the 
competitiveness of the SDP method. 



I. INTRODUCTION 

Prominent simulation techniques for condensed matter 
systems are sampling algorithms like Quantum Monte 
Carlo (QMC) [1-4] and variational algorithms like the 
density-matrix renormalization- group (DMRG) HI, @, 
other tensor network state (TNS) app roaches [7Hl(j|. or 
variational Monte Carlo (VMC) |ll| - [l4| . For a num- 
ber of interesting classes of systems, like frustrated or 
fermionic systems in D > 1 spatial dimensions, the pow- 
erful QMC technique is hampered by the sign problem 
Hence, such systems are then often studied with 
variational techniques, minimizing the energy within a 
certain class of states as, e.g., TNS of a certain struc- 
ture. The energy expectation value of the obtained state 
is necessarily an upper bound to the exact ground state 
energy. 

In this article, a complementary approach is presented. 
By relaxations of the ground state pro blem, one ob- 
tains semidefinite programs (SDP) 01l9|, which 

can be 

solved efficiently on classical computers and yield lower 
bounds to the ground state energy and corresponding 
approximations to few-particle Green's functions. The 
obtained Green's functions allow, e.g., for the study of 
phase diagrams. As the presented SDP method works 
irrespective of the particle statistics, it provides in par- 
ticular a novel route for the study of strongly correlated 
fermionic and frustrated spin systems for D > 1. The 
method can also be used to judge the quality of varia- 
tional algorithms in situations where exact or QMC re- 
sults are not available for comparison. This is especially 
important for VMC methods [HI, [3 and the recently de- 
veloped variational TNS techniques for fermions in D > 1 

The idea is to specify the system by its fc-point Green's 
functions Q^ k \ For systems of fermions, bosons, or hard- 
core bosons (being equivalent to spins-1/2) in a normal- 
ized state p, they read 
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where {a i , a\}i denote the ladder operators with respect 



to some single-particle basis states [25[. The energy 
expectation value E = Tr pH is a function of the (few- 
particle) Green's functions, i.e., E = E(Q). The exact 
ground state energy would be obtained by minimizing it 
among all representable Green's functions 0, i.e., those 
for which a density operator p exists, such that Eq. (QJ is 
obeyed. To determine whether a given Green's function 
G is representable turns out to be a computationally hard 
problem, the famous N -representability problem p6l - [29| . 
which is QMA-complete [30(. An efficient minimization 
of E is however possible if we relax the constraints on 
the Green's function, then, yielding not the exact ground 
state energy but a lower bound, and not the exact ground 
state Green's function but an approximation. Note that 
minimizing E{G) without any constraints is doomed to 
fail as the energy is linear in G. Manageable constraints 
on G can be constructed by imposing the positivity of the 
expectation values of certain positive-definite observables 
like, e.g., the particle density operator a\a i , and so on. 
Such an approach for fermionic G^ and G^ has been 
successfully applied in quantum chemistry; see, e.g., HH 
and references therein. 

Here, we present a systematic method for the con- 
struction and solution of relaxed ground state problems 
for condensed matter systems. Similar to the approach 
known from quantum chemistry, our constraints enforce 
positive expectation values for operators of the form C'C 
with respect to the Green's functions G. More generally 
than before, we (a) choose the constraint operators C to 
be arbitrary polynomials of degree < K in the ladder op- 
erators, which (b) act on suitably chosen subsets of the 
lattice. We (c) exploit the translation invariance of the 
condensed matter systems to (d) work effectively with 
an infinite number of degrees of freedom and describe 
the systems directly in the thermodynamic limit. The 
method yields coupled constraints for G^\ . . . , G^ 2K ^ 
such that the energy optimization problem attains the 
form of an SDP. On basis of the bipolar theorem [32| . 
we further elucidate the mathematical background of the 
method and further possible reductions to the number of 
constraints. 

An alternative method for calculating lower bounds to 
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the ground state energy is due to Anderson. Anderson 
bounds [33l . |34| are obtained by splitting the Hamiltonian 
H into a sum of subsystem Hamiltonians H m that are 
accessible by exact diagonalization. The Anderson bound 
for the ground state energy E° is then given by the sum 
of the ground state energies E^ of the H m , 

H = E m H m E°>J2 m E? n . (2) 

The computation cost for this bound scales exponentially 
in the spatial supports of the operators H m ; see Appx. [B] 
In contrast, the computation cost for the SDP method 
scales only polynomially in the support of the constraint 
operators. In the prominent systems that we studied with 
moderate computer resources, the SDP method typically 
outperforms the Anderson bound. The SDP method has 
the additional advantage of giving access to translation- 
invariant Green's functions, which can be used to study 
phase diagrams etc.. Furthermore, it is possible to choose 
different spatial supports of the constraint operators C 
depending on their degree in the ladder operators. The 
optimal choice for those spatial supports depends on the 
position in the phase diagram. 



where V denotes the set of all density operators and 1Z 
denotes the set of all representable Green's functions, 

K : = {Q\3p e V : Q a = Trpa V &eA } 

= {Q I QM = 1, &t - Qt, Q[B] > V^ }. (7) 

This equality follows from the fact that the only con- 
straints on a valid density operator p are Tr p = 1, p = p', 
and its positivity p >; which is equivalent to requiring 
Tr pB > for all positive semidefmite operators B >z 0. 

Variational methods proceed from Eq. ([6]) by choosing 
some accessible subset of V. Each variational state from 
such a subset yields an upper bound to the ground state 
energy E° . In contrast, for the SDP method, described in 
the following, one chooses an accessible superset J- of the 
set TZ of representable Green's functions, i.e., relaxes the 
constraints. Minimizing the energy in such a superset 
yields a lower bound to E°. A decisive feature of the 
SDP method is that the minimum energies for the chosen 
supersets T can be found certifiably; sec Appx. [X] 

III. SDP METHOD 



II. GROUNDSTATE PROBLEM 

The following description applies to lattice systems 
of fermions, bosons, and hardcore bosons - each corre- 
sponding to a certain algebra for the ladder operators 
{cii}i [251 ] . Spins-1/2 can be treated by mapping them to 
hardcore bosons via the identification a i = S~ , a[ = Sf , 
and Sf = a\a i — |. The generalization to higher spins is 
straightforward. For each subset f2 of the single-particle 
modes {\i}}, let Aq denote the operator basis of normal- 
ordered monomials of degree k in the ladder operators 
for subsystem O. 



= {c 



i\,\Q<m<k,u&Q] (3) 



Each density operator p corresponds to a representable 
Green's function Q, with its fc-point component given by 



Qf :=Trpa 



for a e Aq. 



(4) 



For the moment, let us choose f2 to be the full system. 
Every linear operator B on the Hilbert space can be ex- 
panded in the basis A := Ufc A k as B = J2aeA Its 
expectation value with respect to a state p is then 



TipB= J2^aB a =: G[B], 



(5) 



where Q is the Green's function of p. In this sense, 
Green's functions are linear functionals on the operators. 
For a Hamiltonian H, the ground state problem reads 



E° = min Tr pH = min Q \H] 

pev gen 



(6) 



Solving the ground state problem in general is 
known to be a computationally hard problem; it is QMA- 
complete 35(. Similarly, determining whether given 
Green's functions are representable, i.e., elements of 1Z 
in Eq. (J7J) , is also a QMA-complete problem [3(| . 

A straightforward way to relax the - apparently too 
demanding - constraints represented by 1Z is to require 
the Green's functional G[B] to be non- negative not for 
all B y 0, but only for operators B of the form B = C^C 
with constraint operators C from some set C. Minimizing 
the energy G[H) with respect to Green's functions G from 
the set 

T c ■= {G I G[Id] = 1, = Gl, GpC] > V e& ,}, 

(8) 

yields a lower bound to the ground state energy 

E° > min G\H\ = min VG}fl f , (9) 

as 7Z C Tc- Imposing more and more constraints by 
enlarging the operator set C, the bound approaches E°, 
and the optimal G approaches the ground state Green's 
function. Determining the optimum in Eq. (J9)) is a 
semidefinite programming problem for the variables G&, 
as G[H] = J2a GaH& is linear in G, and the constraints 
G[&C] > VC e C can be written in the form 



to, 



(10) 



where the Hcrmitian matrices Mg- are completely deter- 
mined by the underlying algebra of the ladder opera- 
tors di and the choice for the constraint operator set 
C. Eqs. © and (TIT))) correspond to a standard form 
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for an SDP (UBS; see also Appx. [£] Eq. (JTOJ) results 
from expanding the constraint operators in the basis A, 
C = *YliftCs<j. This yields the constraints in the form 
C;,G[{aya"}C a » > V C - Expanding the oper- 
ators (&')'&" in the operator basis A (by bringing them 
into normal-ordered form) and using G[&] = Ga yields 
the matrices \Ma\a 1 ,a n and the constraints 

J2 GtC%, [M a ]a>,a>>Ca>> = ^ GaCHlaC > V c 

oo' o" a 

from which Eq. (fT0|) follows. 

IV. THERMODYNAMIC LIMIT AND 
CONSTRAINT OPERATORS 

Let us now turn to the specific case of condensed 
matter systems in the thermodynamic limit. Let the 
Hamiltonian be translation-invariant H = ^2 r T r (h) with 

finite-range interaction terms h and the lattice transla- 
tion operator %. We denote the spatial support of the 
K-point terms in h by A K . For a particle-number conserv- 
ing D-dimcnsional model with 2-body nearest neighbor 
interactions, this means for example | A.2 1 = | A.4 1 = D + l, 
and A K = V K ^2,4- Due to the translation invariance of 
H , we can restrict ourselves to translation-invariant den- 
sity matrices and Green's functions. A constraint opera- 
tor set can be constructed by choosing, for each operator 
order k, a subsystem Qfc of the full lattice such that 

fife/ C tt k Wk,k'>k and A K C V K . (11) 

Every such choice of subsystems and the corresponding 
set of constraint operators 

Cn:=spanAn with Aa~[jA^ k (12) 

k 

defines with Eq. flSJ a set Tc n D 1Z of Green's functions. 
The number of Green's function elements Ga occurring 
as degrees of freedom in the SDP is then given by the 
size I .An I of the operator basis. It depends on the parti- 
cle statistics: For a given subsystem 17 of a bosonic sys- 
tem, \Aq\ grows exponentially in the operator order k. 
For fermions, hardcore bosons, and spins- 1/2, \Aq\ ini- 
tially grows exponentially but drops to zero for k > 2\tt\. 
We always choose some K so that ilk = Vk>K- En- 
larging the subsystems fife, systematically improves the 
solution of Eq. ((9]) and increases the computation cost 
polynomially. For given model and computer resources, 
the optimal choice for the subsystems ilk depends on the 
position in the phase diagram. 

V. SYMMETRIES 

Hamiltonian symmetries, like translation or rotation 
invariance, imply that several Green's function elements 
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Figure 1: (color online) Lower bounds E to the ground state 
energy and approximations to correlators for the XXZ chain 
(|13p . The Anderson bound ([2j was calculated with clusters 
of 25 sites. The subsystems ilk for the SDP method ([9| are 
clusters with sizes |fii|, |fia |, and ISI3 1 as specified in the leg- 
end (K — 3). The Bethe ansatz yields the exact ground state 
energy E° and short-range correlators [3^, H3| • 

Ga can be chosen to be identical (e.g., G~ -t = G~ _t 

dajOj, a O a y-^ 

V xy ) and it is sufficient to use in the SDP only one 
representative for each of the corresponding equivalence 
classes. Further, several Green's function elements can 
be zero (e.g., Ga^a = for particle- number conserving 
models). A corresponding block structure in GaMa 
can be exploited to further reduce the computation cost. 

VI. EXEMPLARY APPLICATIONS 
A. XXZ chain 

Let us first address the spin- 1/2 XXZ model 

h = )_\s/.s-; + §y§» + j z s*s*) (13) 

in one spatial dimension (ID). For J z < —1, the model 
is in a gapped ferromagnetic phase with a fully polar- 
ized ground state. In the region — 1 < J z < 1 there 
is the gapless "XY" phase. For J z > 1, the system is 
in the gapped antiferromagnetic Neel phase. The phase 
transition at J z = —1 is of second order and the one 
at J z = 1 is of Bcrczinsky-Kosterlitz-Thouless type [3|| . 
Using comparable (moderate) computer resources, the 
energy bound obtained from the SDP method © im- 
proves substantially on the Anderson bound ^ in large 
parts of the phase diagram; Fig. [1] For J z = 0, where 
the model corresponds to a system of free fermions, the 
SDP bounds reproduce the exact ground state energy 
to high precision. Employing higher-order Green's func- 
tions tends to improve bounds at larger J z . The obtained 
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Figure 2: (color online) Lower energy bounds and correla- 
tors for the 2D XXZ model, complemented by QMC data 
calculated for a square lattice of 16 x 16 sites with periodic 
boundary conditions, and inverse temperature /3 = 96. The 
QMC energies for J z = 0, 1 coincide with earlier results (♦) 
from Ref. [3^,[4(J. The QMC error bars would be smaller than 
the line width, and another simulation with a 32x32 lattice 
produced visually indistinguishable results. We chose SDP 
constraint subsystems f2i and SI2 of square shape with sizes 
as specified in the legend. 



short-range correlators coincide very well with the exact 
Bcthc Ansatz results [IHEzl. 



B. 2D XXZ model 



Let us now consider the spin- 1/2 XXZ model (I13|) on a 



Anderson 5x5 





square lattice. With the identification a 



s; 



a] a, 



i, it maps to a model of interacting hard- 



core bosons obeying the algebra c^a} — (— l) dii dja, = Si 



As displayed in Fig. [2j the SDP method can yield better 
lower bounds to the ground state energy than the An- 
derson bound. However, with the employed computer 
resources, the improvement is arguably not as impres- 
sive as for the other two studied models. As there is 
no exact solution available, we also simulated the model 
with QMC based on the stochastic series expansion Q, 
showing again that the SDP method gives access to the 
correlation functions. 



C. 2D Hubbard model of spinless fermions 

Finally, Fig. [3] shows how the SDP method improves 
upon Anderson bounds for the 2D Hubbard model of 
spinless fermions 

H = ~\ + h.c.) + U - \){h 3 I) (14) 

W> (id) 



free fermions, exact E » 
SDP [9x9,3x3] 
SDP [16x16,4x4] 
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Figure 3: (color online) Lower energy bounds for the 2D Hub- 
bard model of spinless fermions on a square lattice (|14|) . 



on a square lattice. In this case, except for U = 0, no 
exact results are available and QMC is inefficient due to 
the sign problem [L5l - tl7| . 



VII. BIPOLAR THEOREM 

Often, one is only interested in the single- and two- 
particle Green's functions. However, in the presented 
SDP approach, we also introduce higher Green's func- 
tions to improve the approximation. On the basis of the 
bipolar theorem, one can understand that such higher 
Green's functions represent slack variables. The set of 
representable (k < A)-point Green's functions 

Kk :={G\3pEV: Q^=Trpa V k < K ,azA*} 

is convex, as Q is linear in p and T> is convex. Giving up 
on the (inessential) normalization of the Green's func- 
tions, the set TZk ■= {aQ \ Q G TZk, ol 6 R+} becomes a 
convex cone. For a given scalar product, {■, •), the bipolar 
theorem [32| states that (TZ* K )* = TZk, where 



K* K := {B\(B,G) >0 y gelZK } 



(15) 



is the polar cone of TZk- With the choice (B,Q) := 

J2kJ2a&A k B ° = ^1^]' tnc P olar is tnc convex 
cone of all positive semidefinite operators from Bk '■= 

spanlJ^L A k . Due to the bipolar theorem, TZk is hence 
characterized by TZ* K as 



TZ K = {Q I &t = Q% y», G[B] >0V S 



ol- (16) 



BeB K ,By 

So, to obtain (or approximate) the {k < if)-point Green's 
functions, one needs to consider only (k < A')-point op- 
erators B y 0. In this sense, higher Green's functions 
are slack variables, which are only employed in order to 
bring the ground state problem into the form of an SDP; 
Eqs. © and flO]). We showed how constraints G[B] > 
can be enforced in the SDP, for the case that B = &C 
with constraint operators C that are polynomials of de- 
gree < A/2. However, there are also subspaccs of opera- 
tors C of polynom degree > A/2 such that G[C^C] can 
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be evaluated with the (k < AT)-point Green's functions. 
They can hence be taken into account without introduc- 
ing higher Green's functions. A particularly simple space 
of such operators for a particle-number conserving sys- 
tem is given by C = J2i c i a %\ ■ ■ ■ a i m +h.c: For every odd 
m, GJC^G] can be evaluated without requiring 

VIII. CONCLUSION 

We have presented a method for calculating lower 
bounds to the ground state energy of condensed matter 
systems and approximations to the ground state Green's 
functions. Based on certain relaxations of the ground 
state problem one obtains efficiently solvable SDPs. The 
method can also be used for systems with particles of 
mixed statistics, higher spins, etc., by employing the cor- 
responding operator algebras. Our generic considerations 
on the SDP method carry over to quantum chemistry 
problems. A decisive advantage in the variant for con- 
densed matter systems is, however, that translation in- 
variance and locality can be exploited to systematically 
balance the number of degrees of freedom and the com- 
putation cost. 
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Appendix A: Semidefinite programming 

For the estimation of lower energy bounds we use 
semidefinite programming (SDP) [l^Qjl. The SDP algo- 
rithm iteratively approaches the solution of two different 
but related numerical problems, one corresponding to the 
lower energy bound we are looking for. Under certain cir- 
cumstances, the limiting point is a guaranteed optimum 
of both optimization problems and provides us with a 
certified lower bound for the energy. 

Let Q and M a with a = 1, . . . , m be symmetric real 
matrices and H £ R m . The first, and so called primal 
problem, is to find 

n*:=maxTrQ T A (Al) 

such that Tr MjX = H a Ver. (A2) 

The second, and so called dual problem, is to find 

m 

d* := min V g a H a (A3) 
gen™ ^— ; ' 

a— 1 

m 

such that ^2 9° M ° -QbO. (A4) 

CT=1 



The dual problem is our energy minimization problem if 
we let g a = G„ be the the Green's function, M a = Ms 
the constraint matrices occurring in the positivity con- 
dition Eq. (jlOp . and H the Hamiltonian in its vectorized 
representation. Let also Q = —Mid, without coefficient, 
because Gm = 1 is fixed, as it corresponds to the norm 
of the state wc arc searching for. 

We always have that p* < d* , because the cone of pos- 
itive semidefinite matrices is self-dual, implying X, Y y 
^ Tr X T Y > and, hence, 

m m 

< Tr ( 9*M a - Q) T X = ^ g a Tr M T a X - Tr Q T X 

<T=1 <T=1 

m 

= Y,9*H a -TrQ T A. (A5) 

cr=l 

The strong duality theorem states that if the dual prob- 
lem is strictly feasible, i.e., 3g : J2a=i 9vM a — Q >- 
and d* > — oo, then p* = d* and, in particular, p* and 
d* attain their supremum and infimum, respectively. In 
fact, the relaxed ground state problems that we describe 
in this article are strictly feasible, as one can always con- 
struct a mixed quantum state (corresponding to a partic- 
ular vector g) that yields positive expectation values for 
all positive semidefinite constraint operators under con- 
sideration. Hence, there is no gap between the limiting 
points of the SDP and the solution is the actual infimum 
of the energy minimization problem. 

During the optimization procedure, the program 
searches for data underlying an improved result for d* 
or p* within, or at the boundary of, convex cones. 
One cone is the set {A^X y 0} the other is {Z\Z = 
T^=i9°M<j h Q}- Different from, e.g., the simplex al- 
gorithm known from linear programming, which iterates 
along the boundary of the set of solutions, we make use of 
interior-point solvers, which show a superior performance 
in the context of SDP. These solvers consider modifica- 
tions of the original optimization problem and advance 
along a path within the convex cone. Given certain as- 
sumptions, using interior point methods, the SDP can be 
solved to any desired numerical precision within polyno- 
mial time. 



Appendix B: Anderson bounds 

Anderson bounds [H, HH are determined by split- 
ting the Hamiltonian H into M terms H m for which 
the ground state energies are accessible, for example, by 
choosing the terms H m such that their supports are small 
enough to allow for an exact diagonalization of each term. 

H = Ei=i H m (BI) 

The computation cost for determining the ground state 
energies of the operators H m scales exponentially in the 
size of their spatial supports. One can exploit the sparse- 
ness of the Hamiltonians and access, with state of the art 
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computer resources, subsystem sizes of up to about 25 
spins-1/2. Larger subsystems are difficult to address, as 
the H m are not translation-invariant. With diagonal rep- 
resentations H m = J2 m nEm\ rn i n )( rn > n \ an d the lowest 
energy eigenstates |m, 0), one has 

H m h £ m ,„ E °J m > n )( m > n \ = E l ' 1 (B2) 

^HhE m E° m -h (B3) 



i.e., J2 m is a lower bound to the ground state energy 
of the Hamiltonian H . 



M. Suzuki, S. Miyashita, and A. Kuroda, Prog. Theor. 
Phys. 58, 1377 (1977). 

N. V. Prokof'ev, B. V. Svistunov, and I. S. Tupitsyn, 
JETP Lett. 64, 911 (1996). 

O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 
046701 (2002). 

F. Alet, S. Wessel, and M. Troyer, Phys. Rev. E 71, 
036706 (2005). 

S. R. White, Phys. Rev. Lett. 69, 2863 (1992). 

U. Schollwock, Rev. Mod. Phys. 77, 259 (2005). 

H. Niggemann, A. Kfiimper, and J. Zittartz, Z. Phys. B 

104, 103 (1997). 

T. Nishino, K. Okunishi, Y. Hieida, N. Maeshima, and 
Y. Akutsu, Nucl. Phys. B 575, 504 (2000). 



[21 



F. Verstraete and J. I. Cirac, arXiv:cond-mat/0407066 

(2004) . 

G. Vidal, Phys. Rev. Lett. 99, 220405 (2007). 
W. L. McMillan, Phys. Rev. 138, A442 (1965). 

D. Ceperley, G. V. Chester, and M. H. Kalos, Phys. Rev. 
B 16, 3081 (1977). 

F. Mezzacapo, N. Schuch, M. Boninsegni, and J. I. Cirac, 
New J. Phys. 11, 083026 (2009). 

H. J. Changlani, J. M. Kinder, C. J. Umrigar, and G. K.- 
L. Chan, Phys. Rev. B 80, 245116 (2009). 

J. E. Hirsch, R. L. Sugar, D. J. Scalapino, and R. 
Blankenbecler, Phys. Rev. B 26, 5033 (1982). 
M. Takasu, S. Miyashita, and M. Suzuki, Prog. Theor. 
Phys. 75, 1254 (1986). 

M. Troyer and U.-J. Wiese, Phys. Rev. Lett. 94, 170201 

(2005) . 

L. Vandenberghe and S. Boyd, SIAM Rev. 38, 49 (1994). 
F. Alizadeh, SIAM J. Optim. 5, 13 (1993). 
C. V. Kraus, N. Schuch, F. Verstraete, and J. I. Cirac, 
Phys. Rev. A 81, 052338 (2010). 

P. Corboz, G. Evenbly, F. Verstraete, and G. Vidal, Phys. 



[22 

[23 

[24 
[25 

[26 
[27 
[28 
[29 

[30 

[31 
[32 

[33 
[34 

[35 

[36 

[37 

[38 



[39 
[40 



Rev. A 81, 010303(R) (2010). 

C. Pineda, T. Barthel, and J. Eisert, Phys. Rev. A 81, 
050303(R) (2010). 

T. Barthel, C. Pineda, and J. Eisert, Phys. Rev. A 80, 
042333 (2009). 

P. Corboz and G. Vidal, Phys. Rev. B 80, 165129 (2009). 
J. W. Negele and H. Orland, Quantum Many-Particle 
Systems (Perseus Books, Reading, 1988). 
R. H. Tredgold, Phys. Rev. 105, 1421 (1957). 

C. A. Coulson, Rev. Mod. Phys. 32, 170 (1960). 
A. J. Coleman, Rev. Mod. Phys. 35, 668 (1963). 

M. Deza and M. Laurent, J. Comput. Appl. Math. 55, 
217 (1994). 

Y.-K. Liu, M. Christandl, and F. Verstraete, Phys. Rev. 
Lett. 98, 110503 (2007). 

D. A. Mazziotti, Acc. Chem. Res. 39, 207 (2006). 

R. T. Rockafellar, Convex Analysis, Princeton Mathe- 
matical Series (Princeton University Press, Princeton, 
1970). 

P. W. Anderson, Phys. Rev. 83, 1260 (1951). 

T. Wittmann and J. Stolze, Phys. Rev. B 48, 3479 

(1993). 

J. Kempe, A. Kitaev, and O. Regev, SIAM J. Comput. 
35, (2006). 

J. D. Cloizeaux and M. Gaudin, J. Math. Phys. 7, 1384 
(1966). 

G. Kato, M. Shiroishi, M. Takahashi, and K. Sakai, J. 

Phys. A: Math. Gen. 36, L337 (2003). 

Quantum magnetism, Lecture notes in physics, edited by 

U. Schollwock, J. Richter, D. J. J. Farnell, and R. F. 

Bishop (Springer, Berlin, 2004). 

A. W. Sandvik, Phys. Rev. B 56, 11678 (1997). 

A. W. Sandvik and C. J. Hamer, Phys. Rev. B 60, 6588 

(1999). 



